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A  boundary  variation  method  for  the  forward  modeling  of  multi-layered 
diffraction  optics  is  presented.  The  approach  enables  for  fast  and  high-order 
accurate  modeling  of  infinite  periodic  and  finite  aperiodic  transmission  optics, 
consisting  of  an  arbitrary  number  of  materials  and  interfaces  of  general  shape, 
subject  to  plane  wave  illumination  or,  by  solving  a  sequence  of  problems,  il¬ 
lumination  by  beams  The  key  elements  of  the  algorithm  are  discussed  as  are 
details  of  an  efficient  implementation.  Numerous  comparisons  with  exact  so¬ 
lutions  and  highly  accurate  direct  solutions  confirms  the  accuracy,  versatility, 
and  efficiency  of  the  proposed  method. 
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1.  Introduction 

The  ability  to  accurately  and  efficiently  model  multi-layered  diffraction  dominated 

optics  has  numerous  applications,  i.e. ,  modeling  of  integrated  transmission  optics, 
-1*  Corresponding  author:  Jan.Hesthaven@Brown.edu 
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design  and  analysis  of  Bragg  mirrors,  and  integrated  cavities.  Such  methods  would  also 
enable  the  modeling  of  random  variations  in  design  due  to  manufacturing  limitations, 
as  well  as  serve  as  fast  forward  solvers  in  inversion  as  part  of  an  quality  control 
process. 

For  binary  periodic  structures,  the  rigorous  coupled-wave  analysis  has  been  widely 
used  to  model  such  optical  elements  during  the  last  decade.9  However,  the  need  to 
analyze  structures  with  aperiodic  non-binary  features  of  finite  extent  and  bean  il¬ 
lumination  requires  that  alternatives  be  sought.  This  has  lead  to  the  development 
of  several  alternative  approaches,  e.g.,  Finite- Element, 14  Boundary  element,12  Finite- 
Difference  Time-Domain,16  and  Spectral  Collocation.8,10,11  The  latter  two  methods 
both  compute  a  direct  solution  of  the  time-domain  vectorial  Maxwell  equations  and 
are  very  general  in  being  adaptable  to  a  wide  variety  of  geometries  and  physical  set¬ 
tings.  As  the  need  to  model  problems  of  realistic  size  and  complexity  becomes  more 
pressing,  however,  the  memory  and  computational  time  requirements  of  such  direct 
volume  methods  quickly  becomes  a  limiting  factor  not  only  for  the  design  process  but 
also  for  the  analyses  of  particular  structures. 

In  this  work  we  take  a  different  approach,  building  on  the  boundary  variation 
method  proposed  in.3,4,5  These  methods  are  founded  on  the  observation  that  solu¬ 
tions  to  electromagnetic  diffraction  by  a  periodic  structure  depend  analytically  on  a 
variation  of  the  interface.  In  other  words,  diffraction  from  a  smooth  grating  can  be 
determined  from  knowledge  of  reflection  and  refraction  at  a  plane  interface.  Using  this 
result,  fast  and  accurate  high-order  perturbation  schemes  for  finite-size  perturbations 
have  been  developed  for  modeling  of  two-  and  three-dimensional  metallic  and  trans- 
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mission  gratings3, 4,5  illuminated  by  plane  waves.  These  methods  were  subsequently 
extended  to  include  problems  illuminated  by  guided  waves6, 7  and  verified  extensively 
by  comparisons  with  direct  high-order  solution  of  Maxwell  equations. 

We  continue  the  development  of  these  methods  by  demonstrating  the  use  of 
the  boundary  variation  methods  in  an  iterative  approach  to  accurately  and  efficiently 
model  multi-layered  optics,  e.g.,  transmission  optics  where  it  is  essential  to  accurately 
account  for  the  internal  reflections.  To  make  this  feasible,  we  discuss  in  some  detail  the 
implementation  which  relies  on  properties  of  the  scattering  process  to  make  it  efficient. 
The  proposed  algorithm  is  shown  to  perform  well  for  period  and  nonperiodic  problems 
and  results  agrees  very  well  with  directly  computed  reference  solutions,  albeit  obtained 
a  dramatically  reduced  computational  cost. 

What  remains  of  this  paper  is  organized  as  follows.  In  Sec.  2  we  discuss  the 
basic  setup  and  the  essential  details  of  the  formulation.  This  sets  the  stage  for  Sec. 
3  where  we  introduce  the  boundary  variation  method,  first  for  a  single  interface  and 
subsequently  for  the  general  multiple  interface  problem.  Key  elements  of  an  efficient 
implementation  is  also  discussed  here.  In  Sec.  4  we  offer  a  number  of  test  cases  to 
illustrate  the  efficiency  and  capability  of  the  proposed  approach,  while  Sec.  5  concludes 
with  a  few  remarks  and  outlines  future  work. 
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2.  Problem  Setup  and  Model 


As  illustrated  in  Fig.  1  we  consider  a  two-dimensional  situation  in  which  a  monochro¬ 
matic  plane  wave 

EiDC(x,t)  Ae 

=  exp  i(fcinc  •  x  —  ut )  , 

Hiuc(x,t)\  [Ah  \ 

propagating  in  the  homogeneous  region,  fl°,  impinges  on  a  structure  of  multi-layered 
non- magnetic  homogeneous  regions  fl1, . . . ,  QN ,  separated  by  the  smooth  periodic 
interfaces,  T1, . . . ,  TN ,  described  by  the  functions,  /1(x), . . . ,  fN(x)  For  each  region 
Q3  we  have  the  associated  permittivities,  £3 .  Furthermore,  we  have  introduced  the 
complex  amplitudes,  AE  and  AH,  for  the  electric  and  magnetic  field,  respectively, 
while  fcinc  signifies  the  wavevector  of  the  incoming  field,  restricted  to  propagate  in  the 
(x,y)- plane,  i.e.,  fcmc  =  (k™c,  k™).  Recall  that 

Xv  Xv 

Without  loss  of  generality  we  can  normalize  length  and  time  such  that  the  vacuum 
wavelength,  =  1,  and  the  vacuum  speed  of  light  is  cv  —  1. 

The  illumination  of  the  structure  generates  the  fields  (EJ ,  HJ)  in  the  region  Q3  for 
all  j  G  {0,  lf . . .  .N}.  Clearly,  these  fields  satisfy  the  time-harmonic  Maxwell  equations 

VxF  =  iuHJ  ,  V  x  Hj  =  - iuoejE 3  ,  Vj  e  {0.1,,..  .N}  , 

under  the  constraint  that  the  fields  are  solenoidal,  i.e., 

V  •  Ej  =  V  •  Hj  =  0  Vj  G  {0, 1, . . .  .N}  . 
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In  0°,  the  total  field  is  given  as  (Emc,Hmc)  +  ( E°,H° ),  i.e.,  in  this  case  (E3,H3) 
represents  the  scattered  field  while  (E3,H3)  signifies  the  total  fields  otherwise. 

The  fields  in  the  homogeneous  regions,  fh7-1,  and  IP  are  connected  through  the 
boundary  conditions  along  T-7  for  all  j  G  {1, . . . ,  N}.  In  particular,  along  the  interface 
E7,  endowed  with  an  outward  pointing  normal  vector,  h3 ,  separating  the  two  di¬ 
electric  regions  and  iV ,  continuity  of  the  tangential  field  components  requires 


E3~l  +  £inc 

=  h3  x 

1 

i _ 

Hjl  +  ^_i,0  Hinc 

i 

■o» 

tcj 

In  the  special  case  where  T-7  is  assumed  to  be  a  perfect  conductor,  this  condition 
degenerates  to  a  Diric-hlet  condition  on  the  electric  field  as 

hj  x  E3  1  =  ~nj  x  rt,  i,,f;mK  . 

When  solving  Maxwell’s  equation  it  is  often  advantageous  to  express  the  boundary 
condition  on  H3  through  the  condition  on  E3  and  Maxwell’s  equations  themselves  as 
a  Neumann  condition 

h3  x  V  x  +  8^1:0Hmc)  =  -hj  ■  V  ( H3~ 1  +  8^lfiHinc) 

dn3  ’ 

along  the  perfectly  conducting  interface  T-7. 

To  simplify  matters  further,  we  restrict  ourselves  to  problems  where  Q3  can  be 
considered  homogeneous  and  where  the  incoming  fields  are  two-dimensional  and 
either  TE  or  TM  polarized.  This  implies  that  the  local  solutions  to  the  scatte¬ 
ring/penetration  problems  satisfies  the  homogeneous  Helmholtz  equation 

j  =  0,...,N  :  A u3  +  \kj\ V  =  0  ,  (1) 
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where  u 3  =  E3Z  in  the  case  of  TE-polarized  illumination  and  u 3  =  H3  for  a  TM- 
polarized  incoming  wave.  Furthermore,  \k3\  =  ‘l'n\fe3  represents  the  magnitude  of  the 
local  wavevector  under  the  normalization  discussed  in  the  above. 

This  system  of  Helmholtz  equations,  Eq.(l),  must  be  solved  subject  to  the  con¬ 
ditions  that 

u3~x(x,f3(x)) -u3(xj3(x))  =  -8j-lfiumc(xJ3(x))  ,  (2) 

and 

w ))  -  CjJLu'ixJ'ix))  (3) 

on3  onJ 

at  a  general  di-clcctric  interface.  The  constant  Cj  takes  the  values 

1  TE  Polarization 

e3~x fe3  TM  Polarization 
In  the  case  where  r-5  represents  a  perfectly  conducting  object  the  conditions  are 
different  for  the  two  polarizations,  i.e.,  we  have  a  Diric-hlet  condition  in  case  of  a 
TE-polarized  wave 

u3~\xj3(x))  =  ~8j-lfiE™(x,f3(x))  , 
while  we  recover  a  Neumann  condition  on  the  tangential  magnetic  held  as 
E-ui-'(xjUx))  =  Hr(x,fi(x))  , 

along  the  metallic  interface  in  case  the  illuminating  wave  is  TM-polarized. 
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To  complete  the  specification  of  the  problem  we  require  that  the  solutions,  u°  and 
uN ,  are  bounded  at  infinity  and  that  solutions  consists  of  purely  outgoing  waves.  The 
means  by  which  we  enforce  these  conditions  are  closely  related  to  the  computational 
approach  chosen  to  solve  the  above  system  of  coupled  Helmholtz  equations. 

3.  Boundary  Variation  Method 

As  we  solve  the  general  multiple  interface  problem  as  a  sequence  of  single  interface 
problems  we  first  discuss  the  single  interface  scheme  in  detail  and  subsequently  con¬ 
sider  the  extension  to  the  general  case. 

A.  Single  interface  scheme 

Let  us  assume  that  we  only  have  a  single  interface,  T,  separating  the  two  homogeneous 
regions,  and  =  fH,  illuminated  by  a  two-dimensional  TE  or  TM  polarized 

wave.  We  take  the  amplitude  of  the  illuminating  TE  or  TM  wave  is  1. 

The  interface,  T,  is  assumed  to  be  A-periodic  along  x ,  and  described  by  f(x). 
Since  the  incident  wave  is  a  plane  wave,  the  fields,  ■u±(a;),  are  endowed  with  a  similar 
periodicity,  i.e., 

u±(a;  +  L,y)  =  exp(ik1^cL)u±(x,y)  . 

Following,15  we  shall  fix  the  notation  by  introducing 

2vr 

A'  =  —  ,an  =  k™  +  riK  ,  o?n  +  (/^)2  =  l^2  , 

where  K  simply  reflects  the  wavenumber  associated  with  the  periodicity,  an  represents 
the  Bragg  condition  for  lattice  refraction/reflection  while  the  last  condition  expresses 
energy  conservation.  In  accordance  with  the  orientation  of  the  problem,  see  Fig.  1,  we 
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shall  use  the  notation  that  fin  <  0  as  the  incoming  wave  propagates  in  the  negative 
^/-direction.  Also,  to  ensure  boundedness  of  waves  at  infinity,  we  have  that  Im(/3±)  <  0 
for  the  evanescent  waves.  Clearly  there  can  only  be  a  finite  number  of  propagating 
modes  since  for  n  being  sufficiently  large,  fi±  becomes  purely  imaginary. 

Away  from  the  interface,  T,  we  now  express  the  solution,  u±(a:,?/),  as  a  Rayleigh 
expansion 


v; 


:(x,y)  = 


OO 

E 


Bt 


exp 


i(anxTPnV) 


Conservation  of  energy  implies  that 


E  ft 


B 


+ 


+  0)  E  fin 


b: 


=  kl 


nsn+  nen- 

where  II^1  represent  the  subset  of  /3±  corresponding  to  the  propagating  waves,  i.e., 
fin  —  0  is  real.  For  purely  metallic  scattering  the  second  part  of  the  sum  drops  out. 

In  this  setting,  the  unknowns  are  the  Rayleigh  coefficients,  B±,  which  depends 
on  the  profile,  f(x,d,6).  Following,3  we  introduce  a  new  profile 


fs(x)  =  5f(x,L)  , 


i.e.,  <5  =  0  corresponds  to  a  flat  horizontal  interface,  while  <5  =  1  represents  the  original 
profile. 

The  heart  of  the  boundary  variation  method  introduced  in3,4  is  the  assumption 
that  the  Rayleigh  coefficients,  Bfi(S),  can  be  expressed  as  a  Taylor  series  in  <5, 


OO 

BUS)  =  X 


k= 0 


1  <PB±(0)  k 
k\  ddk 


E  dKnfik 


k= 0 
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i.e.,  we  assume  that  Bfi1  and  hence  u^,  are  analytic  in  the  boundary  variation,  <5. 


The  validity  of  this  is  by  no  means  obvious  but  has  nevertheless  been  established 
rigorously  in.2 

To  sketch  the  way  by  which  one  can  obtain  c4,n,  let  us  consider  the  simplest  case 
in  which  is  a  perfectly  conducting  metallic  object  illuminated  by  a  TE-polarized 
plane  wave,  i.e.,  the  boundary  condition  is 

u+(x,  fs(x))  =  —  exp  ik™cx  +  iklyC6f(x) 

^Frorn  the  Rayleigh  expansion  itself  we  have 

00 

=  E  4,nexp  yi{anx  .  (4) 

<5=0  n=—  00 

We  can,  however,  also  evaluate  the  variation  of  u+  with  respect  to  8  at  y  —  0  using 
the  boundary  condition  as 

=  -(«4r£exP(Sy>v)  (5) 

fr~k  dr~k  r  1  dku+' 

ho  ( r  ~  k) !  dVr~k  W  dSk  J  ^,5=0 

Let  us  introduce  the  Fourier  expansion  of  the  periodic  interface,  f(x,L),  as  well  as 
powers  of  it  as 

jr  crfexp(iKix)  .  (6) 

r-  l=-rF 

Using  Eq.(4)  to  evaluate  the  //-derivative  of  l/k\dku+ /dSk  we  can,  by  combining 
Eqs.(4)-(6),  recover  an  explicit  forward  recurrence  for  the  unknown  expansion  coeffi¬ 
cients,  4,n,  °f  the  Rayleigh  coefficients  on  the  form 


1  dku+ 
k}.  d5k 


dk,n=  -(4f)fcCfc,n  + 


k—  1  min[rF,n+(fc— r)F] 

E  E  Ck_r,n_q(-i(lq)k-rdr,q  . 

r= 0  q=max[— rF,n—  (k— r)F] 
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Hence,  given  all  expansion  coefficients,  dS)t,  and  CSjt,  s  <  k,  —sF  <  t  <  sF ,  we  can 
recover  all  expansion  coefficients  for  s  =  k  by  forward  recurrence. 

Albeit  of  a  more  complicated  form,  similar  recurrences  can  be  derived  for  TM- 
polarized  illumination3  as  well  as  for  TE  and  TM  polarized  illumination  of  a  general 
di-electric  interface.4  Recurrences  for  illumination  by  guided  waves  can  be  found  in.6, ' 

While  this  yields  an  approach  for  computing  the  Taylor  expansion  of  the  Rayleigh 
coefficients 

Bn(5)  =  J2dk,nSk  , 

k= 0 


it  is  generally  not  an  easy  matter  to  evaluate  this  expansion  outside  its  circle  of 
convergence.  To  partially  overcome  this  we  express  the  Taylor  series  by  its  Pade- 
approximation  as 


Bn{6) 


cl o  -t-  cl\5  T  ...  T  O-XjS l 

1  +  h6  +  ...  +  bMSM  ’ 


where  the  coefficients  are  found  by  standard  means  by  requiring  that  it  be  equivalent 
to  the  Taylor  series  for  D  =  M  +  L  +  1.  In  subsequent  discussions  we  shall  take 
M  =  L  =  D/2.  As  it  is  well  known,  Pade-approximations  have  remarkably  better 
convergence  properties  than  Taylor  series1  and  suffices  for  the  present  investigation. 

With  this  in  place,  we  can  compute  the  Rayleigh  coefficients  of  the  global  solution, 
^(Xji/),  at  different  angles  of  incidence,  controlled  by  kmc,  using  only  knowledge  of 
the  Fourier  representation  of  the  interface,  T. 


B.  Multiple  interface  scheme 

Exploiting  the  linearity  of  the  Helmholtz  equation  and  considering  a  geometrical 
optics  series,  the  single  interface  scheme  will  be  used  repeatedly  to  form  the  scheme 
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for  multiple  interfaces.  The  approximate  solution  for  the  single  interface  case  is  given 
as  a  plane  wave  expansion 

<2 

Dn  exP 

n=—p 

where  n  G  {—  p, is  the  set  of  waves  for  which  /3±  <  0  is  real.  The  non- 
propagating  waves  are  disregarded  in  the  single  interface  case  as  we  focus  on  modeling 
of  the  far  field.  For  other  applications  requiring  detailed  near-field  information,  these 
evanescent  waves  could  be  included. 

In  the  multiple  interface  scheme  the  solution  for  the  first  interface,  T1,  is  computed 
first.  This  is  a  sum  of  plane  waves,  Eq.(7),  propagating  away  from  the  interface.  The 
waves  traveling  away  from  the  multi-layer  structure  are  not  propagated  any  further. 
However,  the  waves  traveling  down  into  the  multi-layer  structure  are  propagated  to 
the  next  interface  by  updating  the  phase  by  the  phase  factor  exp (i(3d).  Here  d  is 
the  vertical  distance  between  the  interfaces.  The  single  interface  boundary  variation 
method  is  subsequently  used  on  each  of  these  propagated  waves  at  the  new  interface. 
As  this  process  continues,  a  set  of  plane  waves  in  each  region  is  computed,  the  sum 
of  which  gives  the  approximate  solution  in  that  region. 

By  conservation  of  energy  there  must  be  a  finite  number  of  directions  in  which  the 
waves  can  travel  within  each  region,  i.e.,  the  multiple  interface  scheme  only  evokes 
the  single  interface  scheme  a  finite  number  of  times.  When  the  multiple  interface 
algorithm  starts,  the  single  interface  scheme  is  called  a  finite  number  of  times  for  a 
wave  of  unit  amplitude  in  all  the  possible  wave  directions  to  build  an  efficiency  table. 
When  the  single  interface  scheme  is  needed  a  lookup  in  the  efficiency  table  based  on 


i(anxT  Pud) 


(7) 
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the  direction  of  the  incident  wave  and  the  region  it  lives  in  produces  the  efficiencies 
and  directions  of  the  refracted  and  reflected  waves.  The  amplitudes  of  the  refracted 
and  reflected  waves  are  found  by  multiplying  the  efficiencies  by  the  amplitude  of 
the  incident  wave.  Using  the  fact  that  there  is  a  finite  number  of  possible  waves  the 
method  can  group  waves  to  limit  the  computations.  The  following  section  will  make 
this  need  apparent  in  order  to  make  the  approach  computationally  efficient. 


C.  Implementation 


Let  us  first  consider  the  example  of  a  double  interface  structure,  Fig.  2.  The  struc¬ 
ture  has  a  top  interface,  T1,  which  gives  rise  to  three  diffractive  plane  waves  when 
illuminated.  The  second  interface,  T2,  is  taken  to  be  planar  for  simplicity.  With  the  in¬ 
cident  plane  wave,  Emc(x}t),  and  the  top  interface,  T1,  the  single  interface  boundary 
variation  method  is  used  to  compute  the  solution 

E±(x,y)=  y  E$  , 

n=— 1 


where 


Et  =  Bl 


exp 


i(anx- VPnV) 


Vn  e  {-1,0,1} 


We  assume  that  the  boundary  variation  method  emits  a  solution  which  is  the  summa¬ 
tion  of  three  plane  waves  for  the  region  0°  and  three  plane  waves  for  the  region  Q1 . 
To  track  the  solutions  the  implementation  employs  a  solution  set,  Sn,  associated  with 
each  region,  Qn.  These  are  initialized  to  the  empty  set  for  n  e  {0, 1,  2}.  The  solutions 
are  added  to  their  respective  sets  S°  —  S°  U  and  S1  =  S1  U  {E,”}^=_1. 

The  plane  waves,  {K}L-  , ,  travel  away  from  the  interface  and  requires  no  further 
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consideration  prior  to  evaluating  the  far  field.  The  remaining  waves  {.E“}2=_i  are 
traveling  towards  T2  and  we  propagate  these  waves  by  simply  updating  the  phase. 
For  the  wave,  exp(i/3Zid)EZi,  we  need  to  solve  the  single  interface  problem  at  T2. 
However,  as  T2  is  a  flat  plane  we  use  the  Fresnel’s  equations  to  compute  the  reflected 
and  refracted  waves  exactly,  Ezto  an(l  E-ioi  respectfully.  The  waves  are  added  to 
the  solution  sets  so  S1  =  S1  U  and  S2  =  S2  U  {EZi$}-  The  wave  EZij0  is 

propagating  in  the  negative  y  direction  away  from  the  optical  element  and  need  not 
be  considered  further.  However,  EZfo  is  traveling  towards  T1  and  eventually  interacts 
with  T1  from  below,  i.e. ,  the  interaction  is  computed  using  the  single  interface  algo¬ 
rithm  with  the  — /1(x)  profile.  This  process  of  collecting  the  solution  waves  from  the 
single  interface  boundary  condition  and  propagating  the  waves  in  H1  to  accounted  to 
multiple  internal  reflections  is  continued  as  long  as  needed. 

For  the  general  multi-layer  problem  with  N  interfaces,  we  initialize  our  method 
with  the  region  0°  and  the  incident  wave  Emc  by  adding  (.Emc,  0°,  Tttl)  to  the  wave 
set  W°.  Here  Tttl  signify  time-to-live  as  a  measure  of  how  the  maximum  number  of 
internal  reflections  a  wave  and  its  children  can  undergo.  All  of  the  solution  sets  are 
initialized  to  the  empty  set,  Sn  =  {0}.  Also  let  Dl  =  X^-=0  d7  where  rf7  is  the  thickness 
of  the  layer  hi-7.  With  this  notation,  we  must  choose  d°  =  dN  =  0. 

To  initialize  the  efficiency  table  we  first  note  that  in  any  given  layer  the  local 
Bragg  condition  is  a  =  k™c  +  X/p=i  \kj\2mj  where  rrij  is  an  integer  for  j  =  1,  2, . . . ,  N. 
Also  for  each  region  H-7  we  have  conservation  of  energy 

a2  +  (32  =  \kJ\2  . 
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Recall  that  we  only  consider  the  waves  for  which  /3  is  real,  i.e.  there  are  a  finite  number 
of  waves  in  each  layer  and  we  can  find  all  possible  directions  a  wave  will  be  traveling 
in  each  layer  using  the  Bragg  conditions. 

The  single  interface  boundary  variation  code  can  be  run  for  incident  waves  of 
unit  magnitude  for  each  of  these  possible  directions  in  each  layer  and  the  resulting 
efficiencies  are  stored  in  a  separate  hash  table  for  each  layer.  The  key  to  the  hash  table 
is  the  direction  of  the  incident  wave.  In  the  case  of  an  incident  wave  illumination  a  flat 
plane  interface  the  computations  are  simplified  by  using  Fresnel’s  equations17  instead 
of  the  single  interface  boundary  variation  code. 

Assume  now  that  we  are  considering  a  wave  that  has  undergone  m  bounces  and  is 
scheduled  undergo  the  next  bounce.  First  note  that  E  =  A  exp  [i(ax  +  C/3y)],  where 
C  —  —  1  if  the  wave  is  traveling  in  the  positive  y  direction  and  C  —  1  if  the  wave 
is  traveling  in  the  negative  y  direction.  For  an  element  w  =  (E,  Ql,  Tttl)  from  the 
set  of  active  waves,  Wn,  this  wave  is  propagated  by  updating  the  phase  through 
multiplication  by  exp {ij3dl)  and  using  the  single  interface  method  at  the  particular 
interface  with  the  initial  wave  E  exp(i/3dl).  This  yields  the  solution 

^{x.y)  =  ]T  E±  , 

n=—p 


where 


En  =  Bn  exP  [i(anX  =F  f3n  y)  \  Vn  e  {- p ,  -p  +  l,...,q}  . 

Each  wave  is  added  to  a  set  corresponding  to  each  layer,  i.e.,  Sl  =  S'zU{exp(— C /3 Dma,x^,l+C^ ) E^} 
and  Sl+C  =  Sl+C  U  {exp {-C (5Dm^l+c^)E~}qn=_p. 

Next  we  determine  if  we  should  add  the  newly  computed  waves  to  the  set  Wm+1 
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3  ^ 


of  active  waves.  If  the  Tttl  —  1  =  0  then  the  waves  time  to  live  is  up  and  it  is  not 
propagated  further,  i.e.,  no  waves  are  added  to  the  set.  Alternatively,  if  a  particular 
amplitude  is  below  a  given  threshold,  we,  we  can  decide  not  to  propagate  this  wave 
any  further. 

Furthermore,  if  l  +  C  G  {1,2, ...  ,N}  then  the  waves  ({E~}qn=_p,  Vt~,  Tttl  —  1)  are 
added  to  the  set  Wm+1.  If  C  ^  1  V  l  ^  0  then  the  waves  ({E+}l=_pi  f2+,  Tttl  —  1)  are 
added  to  the  set  Wm+1.  When  waves  are  added  to  the  set  Wm+1  it  is  first  checked 
if  a  wave  with  the  same  direction  and  region  is  already  in  Wm+1.  If  so  then  the 
amplitude  of  the  wave  is  added  to  the  amplitude  of  the  wave  already  in  the  set.  If 
there  is  no  wave  with  the  same  direction  in  the  same  region  then  the  wave  is  simply 
added  to  Wm+1 .  Due  to  the  conservation  of  energy,  there  are  only  a  maximum  of 
M  waves  independent  of  m  in  Wm+1.  Without  this  accounting  for  multiple  waves  at 
similar  direction  of  propagation,  the  number  of  active  waves  would  grow  exponentially, 
rendering  the  computational  work  prohibitive. 

The  end  of  going  from  bounce  m  to  m  +  1  is  reached  when  all  the  waves  are 
in  Wm.  To  get  the  approximate  multi-layer  solution  in  the  region  f ll  we  sum  up  the 
elements  in  Sl. 

It  is  possible  to  use  the  L2  difference  in  solutions  of  different  Tttl,s  at  a  few  given 
^/-values  near  the  interfaces  to  evaluate  the  convergence  as  a  measure  of  accuracy.  In 
each  region,  fP,  there  is  a  solution  set  S'-7  which  holds  a  finite  number  of  waves,  e.g., 
s-7.  We  consider  the  error,  e,  in  energy  by  the  relation  for  the  scattering  efficiencies 


1=1  c  1=1 
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where  ej  =  /3i\Bj\2  /  /3mc  is  the  efficiency  and  Bj  exp  i(ajx  +  Cj fify)  is  a  wave  in  Sf 
Note  that  S°  does  not  include  the  incident  wave  and  Cj  =  —  1, 1  denotes  the  direction 
in  which  the  wave  is  traveling. 

It  is  worth  while  emphasizing  that  a  small  energy  deficiency,  e,  does  not  imply 
convergence  to  the  correct  solution  but  simply  checks  self  consistency. 

4.  Verification  and  Convergence 

In  the  following  we  shall  attempt  to  validate  the  accuracy  and  general  performance 
of  the  proposed  scheme  for  several  different  multi-layer  two-dimensional  optical  el¬ 
ements.  When  available,  we  shall  compare  with  exact  solutions.  However,  for  more 
interesting  cases,  these  are  not  available  and  we  compare  against  independently  ver¬ 
ified  highly  accurate  solutions  obtained  with  a  spectral  multi  domain  time-domain 
code.8,10 

A.  Plane  interfaces 

In  the  case  where  all  interfaces  are  planar  the  solution  is  known  on  analytic  form1 7  as 
vP  =  A\  exp  i(a\x  —  (3{y)  +  A2  exp  i(a2x  +  /32y)  , 

within  each  layer,  assuming  that  the  stack  of  layers  is  illuminated  by  a  plane  wave, 
Hexp  [i{ax  —  /%)].  The  unknown  amplitudes  are  found  by  connecting  the  solutions 
through  the  boundary  conditions  Eqs.(2,3)  with  A®  =  A2  =  0.  For  the  wave  directions 
we  have  a{  =  a2  =  a ,  |a|  +  \(3{\  =  |£d+l|  and  |a|  +  \/32\  =  \kf.  This  forms  a  linear 
system  which  can  be  solved  producing  Aj  for  j  —  0, . . . ,  N  and  l  —  1,2.  We  can  then 
compare  the  amplitude  of  the  solution  given  by  the  exact  solution,  |H7|,  with  the 
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amplitude  of  the  solution  given  by  the  multi-layer  boundary  variation  scheme,  \Bf\, 
for  j  =  0, . . . ,  N  and  l  =  1,  2  as  a  function  of  bounces. 

1.  Thin  plate 

For  the  first  test  we  consider  a  normal  TE  wave  of  unit  amplitude  illuminating  a  thin 
plate  (see  Fig.  3).  We  use  Fresnel’s  equations  to  solve  the  single  interface  problems 
and  take  the  wave  amplitude  threshold  of  we  =  10“15,  i.e.,  all  waves  are  allowed  to 
live  for  a  maximum  number  of  bounces. 

In  Fig.  3  we  see  that  the  solution,  or  rather  the  amplitudes,  converges  expo¬ 
nentially  fast  in  the  number  of  bounces  to  the  exact  solution.  The  importance  of 
accounting  for  the  internal  reflections  is  evident  even  for  this  simple  test  case. 

When  comparing  the  solutions  we  look  at  the  error  in  complex  amplitude  of  the 
waves  in  the  exact  solution  Aj  with  the  computed  solution  Bj.  As  we  increase  the 
number  of  bounces,  i.e.,  Tttl  beyond  22  the  error  remains  constant.  This  is  because  the 
wave  amplitude  threshold  we  is  now  fully  responsible  for  terminative  the  the  waves. 

2.  Thin  plate  stack 

A  further  test  of  accuracy  is  illustrated  in  Fig.  4,  reflecting  a  normal  TE  wave  of  unit 
amplitude  illuminating  four  thin  plates  stacked  on  top  of  each  other.  Again  we  use 
Fresnel’s  equations  to  solve  the  single  interface  problems  since  all  of  the  interfaces  are 
planar  and  choose  the  wave  threshold  to  be  we  =  10-15. 

This  problem  provides  a  test  of  the  collecting  of  waves  being  propagated  into  the 
set  Wm+1  into  a  finite  number  of  waves  independent  of  number  of  bounces  m.  The 
work  to  compute  the  solution  increases  linearly  with  the  number  of  bounces  because 


IT 


of  this  collection  process. 

As  is  clear  from  Fig.  4  there  is  exponentially  fast  convergence  to  the  exact  solution 
as  the  number  of  internal  reflections  increases.  We  can  compute  down  to  machine 
precision  with  90  bounces  in  less  than  two  tenths  of  a  second  on  an  average  desktop 
computer.  For  bounces  greater  than  90  the  wave  tolerance  we  determines  when  the 
waves  stop  propagating 

B.  Single  curved  interfaces 

In  the  following  we  provide  a  few  tests  in  which  one  interface  is  curved  while  all  other 
remain  planar.  For  simplicity,  and  supported  by  the  general  verifications  in  the  above, 
we  just  consider  the  two-interface  case.  All  results  are  compared  to  fully  converged 
results  obtained  with  a  high-order  accurate  full  field  solver.8,10 

1.  Shallow  single  curve 

This  problem  consists  of  a  normal  TE  wave  of  unit  amplitude  illuminating  a  thin 
plate  with  one  sinusoidal  interface  and  a  flat  bottom  as  illustrated  in  Fig.  5.  The 
single  interface  interactions  with  the  sinusoidal  interface,  T1,  are  calculated  using  a 
[31/31]  Pade  approximant.  The  interactions  with  T2  are  calculated  using  Fresnel’s 
equations. 

The  wavelength-to-period  ratio  of  T1  is  small,  0.16,  providing  a  simple  test  of  the 
multi-layer  boundary  variation  code.  While  building  the  efficiency  table,  the  single 
interface  boundary  variation  scheme  is  used  8  times  to  compute  possible  wave  interac¬ 
tions  with  T1  and  the  maximum  error  in  the  energy  efficiency  for  the  single  interface 
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cases  is  1.33  x  1CT15. 

In  Table  1  we  increased  the  number  of  allowed  internal  bounces  sufficiently  to 
ensure  that  the  wave  tolerance,  we,  is  the  factor  that  controls  when  the  waves  stop 
propagating.  This  table  shows  that  the  energy  conservation  from  the  single  interface 
scheme  is  kept  as  the  number  of  bounces  increase,  i.e.,  there  does  not  appear  to  be 
severe  problems  of  error  accumulation. 

To  check  the  overall  accuracy  of  the  scheme,  we  illustrate  in  Fig.  6  and  Fig.  7 
the  solutions  obtained  from  the  multi-layer  code  with  results  in  a  different  manner. 
We  see  in  Fig.  6  the  importance  of  considering  the  multiple  internal  reflections.  In 
Fig.  7  we  have  increased  the  number  of  internal  reflections  to  ensure  that  waves  are 
propagated  until  their  amplitude  is  below  we  =  10“15,  yielding  excellent  agreement 
with  the  reference  solution. 

As  one  could  expect,  the  multi  layered  boundary  variation  approximation  be¬ 
comes  more  accurate  as  the  distance  from  the  structure  increases.  This  is  a  conse¬ 
quence  of  not  including  the  evanescent  waves  which  remain  in  the  reference  solution. 
The  excellent  agreement  of  the  fields  at  different  values  of  y  confirms  that  the  phases 
of  the  fields  are  correct  to  a  similar  level  of  accuracy.  This  is  confirmed  by  direct 
comparisons  which  we  omit  for  briefness. 

2.  Deeper  single  curve 

To  further  test  the  algorithm,  we  consider  a  the  case  of  a  TE  polarized  wave  illumi¬ 
nating  a  two-interface  problem  in  which  the  top  interface  is  varies  considerably  more 
than  in  the  above  case,  see  Fig.  8.  In  this  case  the  wavelength-to-period  ratio  of  T1 
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is  0.5.  The  single  interface  interactions  with  T1  are  calculated  using  a  [49/49]  Pade 
approximant,  leading  to  a  maximum  error  in  the  efficiency  for  the  single  interface 
cases  of  0(  10-'). 

As  the  number  of  internal  reflections  increased  to  ensure  that  it  does  not  restrict 
the  propagation  of  the  waves,  the  energy  conservation  is  (P(10-6).  Thus,  as  in  the 
above  simpler  case,  we  find  that  the  energy  conservation  of  the  multi-layer  boundary 
variation  scheme  is  largely  controlled  on  the  maximum  error  in  the  efficiency  for  the 
single  interface  cases. 

In  Fig.  9  we  have  increased  the  bounces  large  enough  so  that  all  waves  are  prop¬ 
agated  until  the  amplitude  is  below  we  =  10-15.  The  agreement  with  the  reference 
solution  is  excellent  and  improves  as  one  moves  away  from  the  structure. 

C.  Lens 

As  an  example  of  a  more  general  nonperiodic  problem,  we  consider  a  Gaussian  lens 
being  illuminated  by  a  normal  TE  wave  of  unit  amplitude.  The  setup  is  illustrated  in 
Fig.  10.  Here  /1(x)  =  0.5exp(— x2)  for  x  G  (—7.5,  7.5]  and  periodically  extended.  The 
periodicity  of  T1  is  set  large  enough  so  that  its  effect  on  the  solution  is  small,  thus 
approximating  a  non-periodic  surface.  For  the  Fourier  transform  of  T1,  needed  by  the 
single  interface  boundary  variation  scheme,  we  use  33  modes,  and  the  solutions  are 
calculated  using  a  [37,  37]  Pade  approximant.  Since  T2  is  flat,  Fresnel’s  equations  are 
used  to  calculate  the  scattering  and  penetration. 

The  multi-layer  boundary  variation  solution  in  Fig.  11  is  for  a  sufficiently  high 
number  of  internal  reflections  to  ensure  convergence,  i.e.,  is  given  for  the  amplitude  is 
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less  than  we  =  10  15.  In  this  case  the  error  in  efficiency  is  0.8.  However,  the  comparison 
with  the  direct  solution  reveals  excellent  agreement,  in  particular  away  from  the  lens. 

D.  Double  curved  interfaces 

As  a  final  test  to  illustrate  the  versatility  of  the  developed  scheme,  we  consider  a 
TE  wave  of  unit  amplitude  illuminating  a  thin  structure  with  two  curved  interfaces, 
illustrated  in  Fig.  12.  In  this  case,  both  T1  and  T2  are  sinusoidal  and  a  [49/49]  Pade 
approximant  is  used  to  calculate  the  single  interface  interactions.  The  number  of 
internal  reflections  is  sufficiently  high  to  ensure  that  all  waves  are  propagated  until 
their  amplitude  is  below  w€  =  10~15 

The  results,  shown  in  Fig.  12,  show  excellent  agreement  between  the  multi-layer 
boundary  variation  scheme  and  the  directly  computed  reference  solution.  The  only 
exception  is  y  =  —0.5.  This  can  be  attributed  to  the  strong  evanescent  waves  produced 
by  the  two  curved  interfaces  and  the  high  contrast.  Here  the  error  in  efficiency  is 
£>(10-13). 

5.  Concluding  Remarks 

The  main  purpose  of  this  work  has  been  to  develop  an  efficient  and  accurate  compu¬ 
tational  approach  to  model  multi  layered  optical  elements,  possibly  with  nonperiodic 
profiles  and  arbitrary  profiles  at  each  interface  separating  layers  of  homogeneous  non¬ 
magnetic  magnetic  materials. 

At  the  heart  of  the  scheme  is  a  very  efficient  boundary  variation  method  for 
accurately  solving  the  problem  of  reflection  and  refraction  by  a  single  interface,  be 
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it  material  or  metallic.  The  methods  exploits  the  linearity  of  Helmholtz  equation 
along  with  the  Bragg  condition  to  significantly  reduce  the  amount  of  space  and  time 
needed  to  compute  the  solutions,  consisting  of  numerous  internally  reflected  waves. 
A  detailed  understanding  of  the  scattering  processes  allows  us  to  collect  the  waves 
into  a  finite  set  of  active  waves,  limiting  the  otherwise  exponentially  growing  set  of 
waves.  Following  a  setup  phase,  tens-of-thousands  of  internally  reflected  waves  can 
subsequently  be  computed  a  very  little  additional  cost. 

In  the  cases  presented  here  the  evanescent  waves  produced  by  the  single  interface 
boundary  variation  scheme  are  disregarded  as  we  have  focused  on  the  solution  away 
from  the  grating  structure.  However,  for  other  applications,  one  could  include  these  in 
a  way  similar  to  that  used  for  propagating  waves.  Given  that  the  kernel  is  the  single 
interface  scheme,  the  limitations  is  set  by  this,  i.e. ,  in  double  precision  arithmetic, 
the  depth-to- wavelength  ratio  of  the  variations  of  the  interfaces  can  not  exceed  0(1). 
Using  extended  precision  can  help  on  this  although  the  cost  typically  increases  also. 

To  check  the  accuracy  of  the  multi-layer  boundary  variation  code  it  has  been  ex¬ 
tensively  validated  through  direct  comparisons  with  high-order  accurate  time-domain 
solutions,  showing  excellent  agreement,  at  a  dramatically  reduced  cost  as  compared 
to  the  direct  solution. 

While  these  results  offer  the  first  step  in  the  development  of  a  general  high-order 
accurate  method  for  the  efficient  modeling  of  multi-layer  diffractive  optics,  a  number 
of  important  issues  remain  open.  Straightforward  extensions  include  illumination  by 
Gaussian  beams,  by  solving  a  sequence  of  problems  subject  to  plane  wave  illumination, 
and  the  use  of  the  threshold  we  to  adaptively  control  the  work  and  requested  accu- 
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racy.  Extensions  to  the  three-dimensional  vectorial  case  is  likewise  straightforward, 
following  the  developments  in.5,7  This  sets  the  stage  for  the  use  of  such  methods  as 
efficient  and  accurate  forward  solvers  in  an  optimal  design  loop,  as  demonstrated  for 
the  simplest  metallic  case  in.13  We  hope  to  pursue  such  developments  in  the  near 
future. 
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Fig.  1  Generic  setup  for  scattering  by  a  problem  with  multiple  interfaces. 

Fig.  2  Specific  example  to  illustrate  scheme  for  a  problem  with  two  interfaces. 

Fig.  3  On  the  left  is  shown  the  problem  setup  for  the  test  with  one  plane  layer  while  the 
right  shows  the  decay  of  error  in  field  amplitudes  as  a  function  of  internal  reflections 
or  bounces. 

Fig.  4  On  the  left  is  shown  the  problem  setup  for  the  test  with  a  stack  of  plane  layers 
while  the  right  shows  the  decay  of  error  in  field  amplitudes  as  a  function  of  internal 
reflections  or  bounces. 

Fig.  5  Problem  specification  for  a  single  layer  with  a  shallow  curved  interface. 

Fig.  6  Ez  computed  at  different  heights,  y,  with  y  —  0  corresponding  to  the  vertical 
position  of  the  shallow  curved  interface.  Illumination  is  TE-polarized  at  normal  in¬ 
cidence.  Results  are  shown  for  a  fixed  number  of  internal  bounces,  and  compared  to 
a  highly  accurate  spectral  solution,  illustrating  the  importance  of  accounting  for  the 
multiple  internal  reflections. 

Fig.  7  Ez  computed  at  different  heights,  y,  with  y  =  0  corresponding  to  the  vertical 
position  of  the  shallow  curved  interface.  Illumination  is  TE-polarized  at  normal  inci¬ 
dence.  Results  are  shown  for  converged  solutions  in  terms  of  internal  reflections  and 
compared  to  a  highly  accurate  spectral  solution. 

Fig.  8  Problem  specification  for  a  single  layer  with  a  deep  curved  interface. 

Fig.  9  Ez  computed  at  different  heights,  y ,  with  y  —  0  corresponding  to  the  vertical 
position  of  the  deep  curved  interface.  Illumination  is  TE-polarized  at  normal  inci- 
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dence.  Results  are  shown  for  converged  solutions  in  terms  of  internal  reflections  and 
compared  to  a  highly  accurate  spectral  solution. 

Fig.  10  Problem  specification  for  a  single  layer  with  an  integrated  lens. 

Fig.  11  Ez  computed  at  different  heights,  y,  with  y  —  0  corresponding  to  the  vertical 
position  of  the  integrated  lens.  Illumination  is  TE-polarized  at  normal  incidence.  Re¬ 
sults  are  shown  for  converged  solutions  in  terms  of  internal  reflections  and  compared 
to  a  highly  accurate  spectral  solution. 

Fig.  12  Problem  specification  for  a  single  layer  with  two  curved  interfaces. 

Fig.  13  Ez  computed  at  different  heights,  y ,  with  y  =  0  corresponding  to  the  vertical 
position  of  the  slowly  varying  interface.  Illumination  is  TE-polarized  at  normal  inci¬ 
dence.  Results  are  shown  for  converged  solutions  in  terms  of  internal  reflections  and 
compared  to  a  highly  accurate  spectral  solution. 
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Table.  1  Convergence  of  scattering  efficiencies  and  relation  to  threshold  value, 


used  in  iterative  approach. 


Fig.  2.  Specific  example  to  illustrate  scheme  for  a  problem  with  two  interfaces. 
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Fig.  3.  On  the  left  is  shown  the  problem  setup  for  the  test  with  one  plane  layer 


while  the  right  shows  the  decay  of  error  in  held  amplitudes  as  a  function  of 


internal  reflections  or  bounces. 
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Thin  plate  stack  setup 
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Thin  plate  stack  error  in  wave  amplitude 


Fig.  4.  On  the  left  is  shown  the  problem  setup  for  the  test  with  a  stack  of 
plane  layers  while  the  right  shows  the  decay  of  error  in  field  amplitudes  as  a 
function  of  internal  reflections  or  bounces. 
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Fig.  5.  Problem  specification  for  a  single  layer  with  a  shallow  curved  interface. 


Ez  amp  Ez  amp 


Ez  amp  y=-3  for  shallow  single  curve 


Ez  amp  y=-1  for  shallow  single  curve 


Fig.  6.  Ez  computed  at  different  heights,  y,  with  y  —  0  corresponding  to  the 
vertical  position  of  the  shallow  curved  interface.  Illumination  is  TE-polarized 
at  normal  incidence.  Results  are  shown  for  a  fixed  number  of  internal  bounces, 
and  compared  to  a  highly  accurate  spectral  solution,  illustrating  the  impor¬ 
tance  of  accounting  for  the  multiple  internal  reflections. 
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Ez  amp  Ez  amp 


Ez  amp  y=-3  for  shallow  single  curve 


Ez  amp  y=-1  for  shallow  single  curve 
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Fig.  7.  Ez  computed  at  different  heights,  y,  with  y  —  0  corresponding  to  the 
vertical  position  of  the  shallow  curved  interface.  Illumination  is  TE-polarized 
at  normal  incidence.  Results  are  shown  for  converged  solutions  in  terms  of 
internal  reflections  and  compared  to  a  highly  accurate  spectral  solution. 
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Fig.  8.  Problem  specification  for  a  single  layer  with  a  deep  curved  interface. 
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Table  1.  Convergence  of  scattering  efficiencies  and  relation  to  threshold  value, 
wf ,  used  in  iterative  approach. 
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Ez  amp  Ez  amp 


Ez  amp  y=-3  for  deeper  single  curve 


Ez  amp  y=-1  for  deeper  single  curve 


Fig.  9.  Ez  computed  at  different  heights,  y,  with  y  —  0  corresponding  to  the 
vertical  position  of  the  deep  curved  interface.  Illumination  is  TE-polarized 
at  normal  incidence.  Results  are  shown  for  converged  solutions  in  terms  of 
internal  reflections  and  compared  to  a  highly  accurate  spectral  solution. 
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Lens  setup 
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Fig.  10.  Problem  specification  for  a  single  layer  with  an  integrated  lens. 
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Ez  amp  Ez  amp 


Ez  amp  y=-2  for  lens  Ez  amp  y=-1  for  lens 


Fig.  11.  Ez  computed  at  different  heights,  y,  with  y  =  0  corresponding  to 
the  vertical  position  of  the  integrated  lens.  Illumination  is  TE-polarized  at 
normal  incidence.  Results  are  shown  for  converged  solutions  in  terms  of  internal 
reflections  and  compared  to  a  highly  accurate  spectral  solution. 
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Double  curved  setup 
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Fig.  12.  Problem  specification  for  a  single  layer  with  two  curved  interfaces. 
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Ez  amp  Ez  amp 


Ez  amp  y=-3  for  double  curved  interfaces 


Ez  amp  y=-2  for  double  curved  interfaces 


Ez  amp  y=-0.5  for  double  curved  interfaces  Ez  amp  y=3  for  double  curved  interfaces 


Fig.  13.  Ez  computed  at  different  heights,  y,  with  y  —  0  corresponding  to  the 
vertical  position  of  the  slowly  varying  interface.  Illumination  is  TE-polarized 
at  normal  incidence.  Results  are  shown  for  converged  solutions  in  terms  of 
internal  reflections  and  compared  to  a  highly  accurate  spectral  solution. 
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